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General-rclativistic models of self-gravitating tori around black holes are constructed 
with a self-consistcnt-field method in compactificd coordinates. The numerical code is 
highly accurate and robust, allowing for the construction of models that exactly fill their 
Roche lobe, when a cusp exists. As a first application, we focus on self-consistent models 
with cusp, having different values of constant specific angular momentum. Scaling all 
results with the mass of the black hole, we find evidence that models with constant 
specific angular momentum that can fill their Roche lobe are still limited by I < 4Mbh 
(as is the case for models constructed in a fixed background metric) even for heavy tori. 

1. Introduction 

Massive relativistic accretion tori around black holes (BHs) can form as transient 
structures in several astrophysical scenarios, including the core-collapse of massive 
star£0 and the merger of neutron star (NS) and NS-BH binaries^SH Recent nu- 
merical simulations have demonstrated that the mass of the torus resulting from bi- 
nary NS-NS or BH-NS mergers, which are candidates for the central engine of short 
gamma-ray bursts (GRBs), can be in the range of - O.O1-O.2M l 6 l 7 l 4 l 5 l In addition, 
studies of the stability of accretion tori have revealed that they can be subject to 
several types of non-axisymmetric instabilities in a number of circumstancea^M 10 ! 11 !, 
which can lead to highly variable and unstable accretion rates. Nonaxisymmetric in- 
stabilities have been recently shown to also exist in fully general-relativistic models 
of massive accretion torP^l 

Another important instability, the axisymmetric runaway instabilitjff^ in models 
that fill their Roche lobe, has been investigated extensively. Numerical simulations 
in a time-varying Schwarzschild background showed that the instability is active 
for tori with constant specific angular momentuirfSl. For the same class of tori, the 
adoption of a time-varying Kerr background resulted in longer timescales for the 
instability to grow, as the rotation of the black hole increased^. Adopting a power- 
law distribution of the specific angular momentum with (even small) positive slope 
stabilizes the torus against the runaway instabilityHESHH The instability was 
shown to be active for constant specific angular momentum tori, even when tak- 
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ing self-gravity into accountJM^H For self-gravitating, non-constant specific angular 
momentum tori, the presence of the instability still depends on the slope of the 
specific angular momentum law, with the range shrinking as the black hole rotation 
increase:^. When self-gravity is taken into account and the torus loses more than 
a few percent of its mass due to mass transfer, it was found that the instability 
is active for different specific angular momentum lawd 22 ! 2 ^ , However, more recent, 
fully general-relativistic simulations of some models of self-gravitating tori with ei- 
ther constant or non-constant specific angular momentum did not show signs of the 
instability on a dynamical timescale^H The parameter space to be investigated for 
the conditions of the onset of the runaway instability is still large and there is a 
need for highly accurate self-consistent models, suitable as numerical initial data 
for simulation codes. 

Analytic models of tori in a fixed spacetime were presented in Refs. [551 HS1 
These can serve as an initial guess for a self-consisted-field approach. Indeed, some 
first self-consistent models in quasi-isotropic coordinates were presented in Ref. 1271 
but with a numerical grid truncated at a finite outer radius. A different numerical 
method was used for constructing self-consistent models of uniformly rotating rings 
around black holeiP^. More recently, a new numerical method was presented for 
constructing initial data of self-consistent models that are suitable for simulations 
using the moving puncture method 29 . Since Roche-lobe filling models can have very 
large outer radii, compared to the radius of the central black hole and, in addition, 
accurate gravitational wave extraction in numerical simulations requires the outer 
numerical boundaries to be placed at even larger distances from the source, the task 
of computing accurate initial data on such different length scales is challenging, but 
can be met with a compactified grid. In particular, a global compactified grid can 
be constructed with a re-definition of the radial coordinate, as was done in Ref. [30] 
for models of rotating neutron stars. This type of compactification is adopted here 
in order to remove the limitations of the method presented in Ref. 1271 

The effect of self-gravity on the properties of a black-hole torus spacetime (and 
on the properties of the torus itself) will, of course, directly depend on the ratio of 
torus mass to black hole mass. Here, we present evidence that even for heavy tori, 
models with constant specific angular momentum, /, can fill their Roche lobe only 
if it is limited by the same relation that holds for models constructed in a fixed 
background metric, I < 4Mbh, where Mbh is the mass of the black hole. 

The paper is organized as follows: In Sec. [2] we summarize the main definitions 
and describe the numerical setup. In Sec. [3] we discuss in detail a representative 
model, while Sec. [4] presents the main properties of a sequence of models with 
constant specific angular momentum. A brief discussion is found in Sec. [5j 

Throughout the paper, we assume a signature of (— , +, +, +) and set c = G = 1. 
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2. Tori in stationary, axisymmetric spacetimes 

The details of our formalism and numerical method can be found in Ref. 1311 Here 
we only summarize the main equations: 

We assume that the spacetime is stationary and axisymmetric, with circular 
velocity field, while the matter is described as a perfect fluid with total energy 
density e, rest mass density p and pressure p. We further assume that the torus has 
a constant specific angular momentum I := —u^/ut, where u a is the 4- velocity. The 
hydrostationary equilibrium equation can be written in the form 

^ = - Vln ^ + l^I< (1) 

where f2 := u^/u* is the angular velocity and has a first integral only under specific 
conditions. For a barotropic equation of state of the form p — p(p) one can define 
the log-enthalpy 

with, dH/dp = (e + p) _1 . Then, the l.h.s. of becomes (e + p)~ x Vp = VH and if 
fi = 0,(1), we can write 

VH = -VW, (3) 
where — VW is the effective gravity, given by the gradient of the effective potential 

W = In(-ut) - / - — —dl + const. (4) 

It follows that the surfaces of constant pressure coincide with the equipotential 
surfaces and at the location of maximum density, VH = VW = (the constant in 
Q can be set by requiring that W vanishes at infinity). 

For the case of constant specific angular momentum the first integral of ([T]) 
becomes 

H + ln(— Ut) = const. (5) 
and for the homcntropic, polytropic equation of state (EOS) 

p=Kp T , (6) 

t=P+f& V (7) 

where K is the polytropic constant and T is the polytropic exponent, one obtains 

H = \n[l + r/(T-l)Kp r - 1 }, (8) 

The density distribution is given algebraically, as an implicit function of Q (for given 
EOS, Iq, inner radius rj n and stationary, axisymmetric spacetime): 



r - i 

KT 



(9) 
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For our self-consistent numerical scheme we use the analytic AJS solutions^ as 
an initial guess. In isotropic coordinates (t, r, 9, <fi), the metric of the Schwarzschild 
spacetime is 



ds 2 = 



1 - M/2 n ' ' rN 1 



dt 2 +[l + ^) {dr 2 + r 2 d9 2 + r 2 sin 2 9d(j) 2 ) 1 (10) 



M 

,11- 4- I I -+- 

1 + M/2r 

where the isotropic radial coordinate r is related to the Schwarzschild radial coor- 
dinate *by r = 1/2 [*— M + V z 2 — 2M 2 . Requiring W(r = oo) = 0, the effective 
potential is 

W = \n{-ut). (11) 

A torus of finite size extends between r = r- ln and r = r out in the equatorial 
plane. At the density maximum, V r p = and the fluid elements there move as free 
test particles with Keplerian angular momentum 



v/M(l + M/2r) D r 2 

(1 - M/2r) 2 (M + M 2 /4r + rf' 2 

has a minimum at the marginally stable circular orbit for test particles, r ms = 
1/2(5 + 2\/6) ~ 4.950, which restricts r max > r ms =► I > l ms = 3.674M. The 
radius of the marginally bound orbit is the location where W — for I — 4M, i.e. 
r m b ~ 2.914. The largest of three roots of I — Ik corresponds to r = r max , while the 
intermediate root, r = r cusp , corresponds to the existence of a cusp, where again 
V r p = 0. For a given EOS and black-hole mass M and for a given value of I, there 
exists a one-parameter family of different finite-size tori with r cusp < rj n < r max , 
which all have the same r max . For a chosen r m in this range, one obtains the density 
distribution from ^ . The parameter space of possible equilibrium configurations is 
further limited by the existence of the marginally bound orbit, r m b. For / < l m \, = 
4M, the effective potential W is always negative in the equatorial plane, while for 
I = i m b- W = at r m b (an inflection point) and the location of the cusp is at 
r in = r mb- F° r I > imb there exists a region in the equatorial plane with W > 0. 
Thus, for any I > Z m b, only tori without a cusp are possible, with the inner radius 
limited by the location of W — 0. 

For self-gravitating tori, we use the metric in quasi-isotropic coordinates 

ds 2 = ~e 2v dt 2 + e 2a (dr 2 + r 2 d9 2 ) + e 2 ^ v \ 2 sin 2 0(# - ujdt) 2 , (13) 

where 7, a and to are metric functions that depend only on the coordinates r and 
9 (see Ref. l32|) . Because 7 and v are divergent at the event horizon, we will use, 
instead, the functions^ B := e 7 and A := e", with boundary conditions on the 
event horizon B = 0, A = and uj = u^, where uj^ is the constant angular velocity 
of the horizon. The horizon can be chosen to be a sphere of constant radius r = h . 
On the axis of rotation, the condition = 7—^ ensures local flatness. 

In the angular direction we use the coordinate \l := cos9 instead of 9, while 
in the radial direction we use a compactified, nondimcnsional coordinate s, defined 
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through 

r:=r e ^— , (14) 

where r e is a suitably chosen radius (here we use r e — r out , so that the outer 
edge of the torus in the equatorial plane corresponds to s = 0.5). Through this 
compactification, the infinite domain r — > [0, +oo) is mapped onto the finite domain 
s — > [0, 1]. The location of the horizon is at s = s , so that h = r e s / (1 — so)- 
The field equations for the three metric functions A, B and uj are of elliptic type 

V 2 A=S a (s,m), (15) 

(^+^!-^<'|) B = s »<->' < 16 > 

where V 2 is the flat-space Laplacian, and S\, Sb and S u are source terms. These 
elliptic-type equations are inverted using appropriate Green's functions, while the 
equation for the metric function a is integrated as an ordinary differential equation. 

While we use nondimensional units of c = G = = 1 in the numerical 

procedure, the solution can then be scaled to any desired black hole mass. Our grid 
is equidistant in (s,fi) with 2001 x 1001 grid points. For a chosen polytropic index 
r, the solution space of AJS models is three-dimensional: the parameter / controls 
the rotation of the torus, r e /h controls its size and K (or p max ) its mass (through 
the variation of the equation of state) . 

Our fixed-point-iteration method is described in detail in Ref. 1311 Here, we note 
that our numerical procedure yields directly a self-consistent model, starting from 
the AJS solution as an initial guess. In order to construct models that exactly 
fill their Roche lobe, we must be able to compute models that violate this limit, 
in order to approach it by bisection. For this purpose, we also construct models 
starting with an AJS solution that overfills the Roche lobe, i.e. with a value of 
Win larger than W CU sp, but truncated at r cusp (similar models were constructed in 
other works, in order to initiate mass accretion onto the black hole). Exploring the 
allowed parameter range, we then find models that exactly fill their Roche lobe, up 
to a single grid point. 

3. A representative model 

A representative solution is constructed for a polytropic equation of state with in- 
dex N = 3 and constant if/Afg^ =17.34, a ratio of outer torus radius to horizon 
radius r e /h =A9 and with = 0. The specific angular momentum in the initial 
guess is Z/Mg^ s =3.675, while it is Z/Mbh=3-911 in the converged, self-consistent 
solution. A detailed list of various properties, scaled by the mass of the black 
hole, is displayed in Table [T] (see Ref. [3T] for definitions) . The ratio of torus mass 
to black hole mass is Mt/Mbh=0-21, while the total mass of the spacetime is 
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Table 1. Properties of a representative self-consistent model of 
a BH-torus system, with N = 3, «/ M BH=3.911, T- out /h =49.005, 
Wh = and if/Aifg^ V =17.339. All radii are coordinate radii in the 
quasi-isotropic metric. 



Asymptotic spacctime mass 


M/M BH 


1.1730 




Torus mass 


Mt/Mbh 


0.2097 




Maximum density 




6.6126x10" 


-5 


Inner torus radius 




3.0413 




Radius of maximum density 




8.1931 




Outer torus radius 


fout/MBH 


23.6400 




Torus rest mass 


M /M BH 


0.1964 




Torus internal energy 


Ut/M B h 


2.9878x10" 


-3 


Torus rotational energy 


Tt/Mbh 


9.5559x10" 


-3 


Torus gravitational potential energy 


Wt/Mbh 


-3.5865x10" 


2 


Torus rotational to potential energy 


\T/W\ T 


0.2664 




Torus angular momentum 


Jt/mIk 


0.7390 




Black hole horizon radius 


h /M BH 


0.4824 




Black hole Komar charge 


Mh/Mbh 


0.9633 





M/Mbh=1-17. ln these coordinates, the horizon of the black hole is at a radius of 
/io/Mbh=0.482 (compared to 0.5 for the AJS solution). The torus itself extends from 
r in /M B H=3.04 to r out /M B H=23.64. The relation M h /M B h < (M - M )/M BH < 1 
holds as expected^. 

Notice that, the gravitational potential energy of the torus is 

W T := M - M BH - M - T T - U T , (18) 

and includes both self-energy and binding energy between the black hole and the 
torus. While all previous quantities are determined with high accuracy (as was 
checked by proper convergence studies with grids of different sizes), this is not the 
case for Wt, since it results as the difference between considerably larger quanti- 
ties. Even for the largest-size grid we can afford at present, a convergence study 
showed that the numerical error in Wt can range from a few % up to roughly 10%, 
depending on the mass of the torus. This affects also the ratio |T/W|t- We will 
thus report these two quantities only up to such an uncertainty. 

The four panels in Fig. [I] show the distribution of the four metric functions 
a, B, A and w, as a function of the nondimensional compactified coordinate s 
in the equatorial plane (we remind here that the outer radius of the torus is at 
s = 0.5). The first three functions (solid lines) are compared to the corresponding 
Schwarzschild solution. The largest impact of the self-gravity of the torus is in 
the metric function A, while B remains largely unaffected. The metric function u> 
vanishes for the Schwarzschild metric, while in the self-consistent solution it shows 
a peak inside the torus. At the horizon and at infinity, all metric functions attain 
their expected boundary values, enforced by the Green's functions approach. 

Fig. [2] shows the isocontours of the effective potential W in Cartesian coordi- 
nates (X, Z) scaled by the outer radius of the torus. For this particular model the 
marginally bound isocontour, W = 0, does not cross the equatorial plane. Instead, 
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s s 

Fig. 1. Distribution of the four metric functions of the representative model of Table [l] in the 
equatorial plane (solid lines), as a function of the nondimensional compactificd coordinate s (the 
outer edge of the torus is at s = 0.5). The horizon of the black hole is at so = 0.02. The corre- 
sponding Schwarzschild solution (the spacetime used in the AJS limit) is shown as dashed lines. 



all isocontours that cross the equatorial plane have W < and a cusp exists. Matter 
fills exactly the Roche lobe (thick black line), which has a cusp in the equatorial 
plane. Inside the Roche lobe, the isocontours of the effective potential W coincide 
with corresponding isodensity contours. 

Fig.[3]shows the distribution of the effective potential W in the equatorial plane, 
as well as the difference of the Keplerian specific angular momentum Zk from the 
constant value I inside the torus. For this model with cusp, the inner edge of the 
torus coincides with a local maximum in W and with the location where — I = 0. 
The close agreement between these two quantities demonstrates the high accuracy 
of the numerical solution. The radius of the marginally stable orbit, r ms , corresponds 
to the local minimum in — /, while the radius where the density has a maximum, 
r max , coincides with a local minimum in W. Furthermore, at the outer radius of the 
torus, r out , the effective potential has the same value as at the cusp. 
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X/fout 



Fig. 2. Isocontours of the effective potential W, in Cartesian coordinates, scaled by the outer 
radius of the torus, r ou t, for the representative model of Table [l] The torus fills exactly its Roche 
lobe (thick black line), which has a cusp in the equatorial plane. Inside the Roche lobe, the 
isocontours of W and of the density, p, coincide (thin black lines). Isocontours of W outside 
the torus can have negative or positive values (dashed lines), separated by the marginally-bound 
surface (W = 0). The horizon of the black hole is a sphere with radius ho/r ou t = 0.02 (not shown 
here) . 



4. A sequence of models with cusp 

In the AJS limit, models with constant specific angular momentum are required 
to have I > 3.674Mg^ s , in order for a finite-size torus to exist. We construct a 
sequence of self-consistent models with different torus mass, that all fill exactly their 
Roche lobe, with I > 3.675M^ S . Tabled summarizes the main properties of models 
belonging to this sequence, with the first line corresponding to the AJS solution (no 
values for Mt and Jx/Mg H are given, since the AJS solution is independent of the 
torus mass). In the AJS limit, the torus has a very small size, extending between 
Ti n /Msu = 4.84 and r out /MBH = 5.24. In the self-consistent solution, the value 
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-0.05 



-0.1 
^ -0.15 
-0.2 



-0.25 



-0.3 

0.2 0.4 0.6 0.8 1 1.2 

r/M BH 

Fig. 3. The effective potential W and — I in the equatorial plane for the representative model 
of Table [I] The radius of the inner edge of the torus, 7*i n , the marginally stable orbit r ms , the 
location of the density maximum, r max , and the radius of the outer edge of the torus, r out , are 
shown with dotted vertical lines. A dotted horizontal line connects the inner and outer edges of 
the torus at the same value of the effective potential. 



Tabic 2. Properties of a sequence of models with a cusp with N = 3, I /M^ =3.675 
and uih = 0. 



Mt/Mbh 


M/M BH 


l/M Bn 


n n /M BH 


fmax/A'/BH 


r ou t/M B H 


Jt/m£ h 




1.000 


3.675 


4.838 


5.092 


5.244 




0.134 


1.107 


3.859 


3.193 


7.763 


18.235 


0.460 


0.296 


1.249 


3.948 


2.888 


8.765 


39.329 


1.079 


0.366 


1.312 


3.952 


2.784 


9.025 


58.745 


1.357 


0.438 


1.379 


3.947 


2.723 


9.167 


85.049 


1.644 


0.546 


1.484 


3.948 


2.660 


9.293 


133.360 


2.074 



of I/Mbh changes, but in the last four models of the sequence it stays practically 
constant at l/M B n = 3.95. Increasing the mass of the torus only increases the total 
mass-energy of the spacetime, but not the mass of the BH. Along the sequence, 
the radius of the inner edge of the torus decreases, approaching the horizon, while 
the radius of the maximum density increases. The most dramatic increase is in the 
outer radius of the torus, which grows more than 25 times for the most massive 
model considered here, compared to the AJS limit. Part of this difference can be 
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Fig. 4. Outer surface for each of the models presented in Table [2] (except for the AJS limit) 
in Cartesian coordinates. With increasing torus mass, the outer edge of the torus, as well as its 
vertical extent, increase significantly. The horizon of the black hole in each case is a sphere with 
radius ranging from ho/MsH = 0.47 (largest torus) to ho/M^n = 0.49 (smallest torus) and is not 
shown here. 

attributed to the fact that I/Mbh increases. Nevertheless, the last four models in 
Table [2] constitute a proper sequence of models with practically the same Z/Mbh, so 
they can be compared on an equal footing. Among these models, the outer radius 
also increases significantly with increasing mass of the torus. This can been seen in 
more detail in Fig. |1J which compares the structure of the different tori in Table [2] 
(showing only the surface for each model), in Cartesian coordinates scaled by Mbh- 
The total vertical extent of the torus approaches twice its extent in the equatorial 
plane, as the torus mass increases. 

Fig. [5] shows the effective potential in the equatorial plane, for the sequence of 
models in Table [2] As the torus mass increases, the value of the effective potential 
at the inner edge of the torus, W[ n , (shown with a dashed line connecting the 
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-0.01 
-0.02 



-0.03 
-0.04 
-0.05 



-0.06 



-0.07 

2 4 6 8 10 12 

r/M BH 

Fig. 5. The effective potential in the equatorial plane, for the sequence of models presented in 
Table [2] The shape of the potential changes significantly with increasing torus mass, due to self- 
gravity effects. A dashed line connects the local maximum for each model, that corresponds to r; n , 
while another dashed line connects the local minimum for each model, that corresponds to r ma x- 

local maxima in the effective potential) increases significantly, which can be seen 
m more detail in Fig. |6j It appears that W- m only asymptotically approaches the 
marginally-bound value of W — and does not cross this limit. This is natural, since 
marginally bound tori close their Roche lobe at infinity. Thus, a sequence of bound 
models with increasing mass would first need to reach a marginally bound state 
with infinite radius, before connecting to a sequence of models that have no cusp. 
We conclude that for tori with constant specific angular momentum, the limitation 
of I < 4Mbh may also hold for self-consistent models, as it does in the AJS limit, 
although a more detailed investigation is needed to locate the precise limit. 

5. Discussion 

We present an improved method for constructing general-relativistic models of self- 
gravitating tori around black holes. It is a self-consistent-field method, with the 
field equations inverted using Green's functions and its improvement over an existing 
approach^ consists in the use of compactified coordinates. This allows for the metric 



i i | i i i | i i i | i i i | i i i | i i r 

AJS limit 
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0.1 0.2 0.3 0.4 0.5 

Mj/Mbh 



Fig. 6. The value of the effective potential Wj n at the inner edge of the torus for the sequence of 
models presented in Table [2] With increasing torus mass, Wi n approaches, but does not cross the 
marginally-bound value of W- ln = 0. 



to be computed up to spatial infinity, resulting in highly accurate initial data for 
numerical evolution codes. These initial data have already been used successfully in 
studying nonaxisymmetric instabilities in fully general-relativistic models of massive 
accretion torP^l. Of interest is also a more detailed study of the precise conditions 
under which the axisymmetric runaway instability is active, since the numerical 
code is highly accurate and robust, allowing for the construction of models that 
exactly fill their Roche lobe, when a cusp exists. 

Scaling all results with the mass of the black hole, we found evidence that models 
with constant specific angular momentum that can fill their Roche lobe are still 
limited by I < 4Mbh ( a s is the case for models constructed in a fixed background 
metric) even for heavy tori. The precise limit must be determined with a more 
detailed study. More generally, we plan to investigate the allowed parameter space 
for different types of self-gravitating tori, including the case of rotating black holes 
and more general rotation laws and equations of state. 
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